Generating final model predictions after relativizing
Alice Stears
2025-11-14
In this document, we integrate our ecoregion classification model as well as models for each different functional group and ecoregion we considered into CONUS-wide predictions of both absolute and relative cover. We show maps of these model predictions that are generated using contemporary climate, weather, and soils data, and assess how these predictions compare to observed data using maps of residuals, quantile plots, bias, and RMSE. We also show model predictions of both absolute and relative cover made using modeled climate and weather data from the end of the 20th century under RCP 8.5. We used data from two different generalized circulation models (GCMs) from CMIP5: “BNU-ESM”, a cooler and wetter model, and “IPSL-CM5A-MR (France)”, a hotter and drier model. For these future predictions of cover, we display maps that show how these future model predictions compare to modeled contemporary cover. We present all of this information (maps of contemporary cover predictions, residuals, quantile plots, RMSE, bias, and maps of future cover based on GCMs) for both absolute and relative cover.
Setup the environment:
library(tidyverse)
library(sf)
library(terra)
library(kableExtra)
library(knitr)
library(USA.state.boundaries)
library(tidyterra)
library(ggpubr)
library(USA.state.boundaries)
library(ggrepel)
# set ggplot theme default
theme_set(theme_classic())
# Load User defined parameters and functions
print(params)## $readParams
## [1] TRUE# set to true if want to run for a limited number of rows (i.e. for code testing)
readParams <- params$readParams
## function to maek predictionss f
makePredictions <- function(predictionDF, modelObject, scale = TRUE) {
##
# predictionDF <- climDatPred
# modelObject <- bestLambdaMod_grassShrub_totalHerb
# # ##
# predictionDF = modDat_quantile,
# modelObject = forb_CONUS, scale = FALSE
#
# pretend to scale the input variables so the model object can predict accurately
if(scale == TRUE) {
predictionDF <- predictionDF %>%
mutate(across(all_of(prednames), base::scale,scale = FALSE, center = FALSE))
}
# modelPredictions
modelPreds <- predict(object = modelObject, newdata = predictionDF, type = "response")
# add predictions back into the input data.frame
predictionDF <- predictionDF %>%
cbind(modelPreds)
# truncate all predictions to max out at 100
#predictionDF[predictionDF$modelPreds>100 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 100
predictionDF[predictionDF$modelPreds < 0 & !is.na(predictionDF$modelPreds),"modelPreds"] <- 0
# print predicted data
return(predictionDF)
}
## source functions for quantile plots
source("../../../Functions/glmTransformsIterates.R")
source("../../../Functions/transformPreds.R")
# Load Data ---------------------------------------------------------------
# data ready for modeling (that has been scaled)
modDat_1_s <- readRDS("./models/scaledModelInputData.rds")
# get the soil raster, which we'll use for rasterizing the imagery
soilRastTemp <- readRDS("../../../Data_processed/SoilsRaster.rds") %>%
terra::unwrap()
# make a map of the predictions
test_rast <- rast("../../../Data_raw/dayMet/rawMonthlyData/orders/70e0da02b9d2d6e8faa8c97d211f3546/Daymet_Monthly_V4R1/data/daymet_v4_prcp_monttl_na_1980.tif") %>%
terra::aggregate(fact = 12, fun = "mean", na.rm = TRUE)
# download map info for visualization
data(state_boundaries_wgs84)
cropped_states <- suppressMessages(state_boundaries_wgs84 %>%
dplyr::filter(NAME!="Hawaii") %>%
dplyr::filter(NAME!="Alaska") %>%
dplyr::filter(NAME!="Puerto Rico") %>%
dplyr::filter(NAME!="American Samoa") %>%
dplyr::filter(NAME!="Guam") %>%
dplyr::filter(NAME!="Commonwealth of the Northern Mariana Islands") %>%
dplyr::filter(NAME!="United States Virgin Islands") %>%
dplyr::filter(NAME!= "U.S. Virgin Islands") %>%
sf::st_sf() %>%
sf::st_transform(sf::st_crs(test_rast)))
#sf::st_crop(sf::st_bbox(modDat_1_sf)+c(-1,-1,1,1))
## add ecoregion boundaries (for our ecoregion level model)
regions <- sf::st_read(dsn = "../../../Data_raw/Level2Ecoregions/", layer = "NA_CEC_Eco_Level2") ## Reading layer `NA_CEC_Eco_Level2' from data source
## `/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/cleanPED/PED_vegClimModels/Data_raw/Level2Ecoregions'
## using driver `ESRI Shapefile'
## Simple feature collection with 2261 features and 8 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4334052 ymin: -3313739 xmax: 3324076 ymax: 4267265
## Projected CRS: Sphere_ARC_INFO_Lambert_Azimuthal_Equal_Arearegions <- regions %>%
st_transform(crs = st_crs(test_rast)) %>%
st_make_valid()
ecoregionLU <- data.frame("NA_L1NAME" = sort(unique(regions$NA_L1NAME)),
"newRegion" = c(NA, "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "dryShrubGrass",
"dryShrubGrass", "Forest", "Forest",
"dryShrubGrass", "Forest", "Forest",
"Forest", "Forest", "dryShrubGrass",
NA
))
goodRegions <- regions %>%
left_join(ecoregionLU)
mapRegions <- goodRegions %>%
filter(!is.na(newRegion)) %>%
group_by(newRegion) %>%
summarise(geometry = sf::st_union(geometry)) %>%
ungroup() %>%
st_simplify(dTolerance = 1000)
## contemporary climate data
climDat_temp <- readRDS("/Users/astears/Documents/Dropbox_static/Work/NAU_USGS_postdoc/PED_vegClimModels/Data_processed/EcoRegion_climSoilData.rds")
climDat_timeTemp <- climDat_temp %>%
mutate(spatialID = paste0(round(x,4), "_", round(y,4))) %>%
drop_na()
spatID_lu <- data.frame("spatialID" = unique(climDat_timeTemp$spatialID),
"uniqueID" = 1:length(unique(climDat_timeTemp$spatialID)))
climDat_timeTemp <- climDat_timeTemp %>%
left_join(spatID_lu)
## because of the spatial joining we did, there are some "locations" that have data for less than the true time span we want (2010-2020), or multiple observations per year... remove those so we have the even representation over the time span we're looking for
set.seed("12011993")
goodTable <- data.frame(table(climDat_timeTemp$uniqueID))
goodIDs_temp <- as.numeric(goodTable[goodTable$Freq == 11,"Var1"])
#goodIDs <- unique(names(table(climDat_timeTemp$uniqueID))[(table(climDat_timeTemp$uniqueID)==11)])
# select 1000 random IDs
goodIDs <- goodIDs_temp[sample(x = 1:length(goodIDs_temp), size = 1000, replace = FALSE)]
## get climate data for 500 different points randomly chosen across CONUS for all years in the contemporary climate/weather/soils dataset
climDat_time <- climDat_timeTemp %>%
filter(uniqueID %in% goodIDs) %>%
dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE,
NA_L1NAME, NA_L1KEY, year, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity, uniqueID, spatialID) %>%
rename("tmin" = tmin_meanAnnAvg_CLIM,
"tmax" = tmax_meanAnnAvg_CLIM, #1
"tmean" = tmean_meanAnnAvg_CLIM,
"prcp" = prcp_meanAnnTotal_CLIM,
"t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
"t_cold" = T_coldestMonth_meanAnnAvg_CLIM,
"prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
"prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM,
"prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
"prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM, #3
"abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM,
"isothermality" = isothermality_meanAnnAvg_CLIM, #4
"annWatDef" = annWaterDeficit_meanAnnAvg_CLIM,
"annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
"VPD_mean" = annVPD_mean_meanAnnAvg_CLIM,
"VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
"VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
"VPD_max_95" = annVPD_max_95percentile_CLIM,
"annWatDef_95" = annWaterDeficit_95percentile_CLIM,
"annWetDegDays_5" = annWetDegDays_5percentile_CLIM,
"frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM,
"frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM,
"soilDepth" = soilDepth, #7
"clay" = surfaceClay_perc,
"sand" = avgSandPerc_acrossDepth, #8
"coarse" = avgCoarsePerc_acrossDepth, #9
"carbon" = avgOrganicCarbonPerc_0_3cm, #10
"AWHC" = totalAvailableWaterHoldingCapacity,
## anomaly variables
tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom, #19
t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom, #22
prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23
prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25
isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom, #28
VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom, #30
VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom, #31
VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33
annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom , #34
frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35
frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
) %>%
dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))
# rename
climDat <- climDat_temp %>%
filter(year == 2016) %>%
dplyr::select(tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_3yrAnom, NA_L1CODE,
NA_L1NAME, NA_L1KEY, newRegion, x, y, soilDepth:totalAvailableWaterHoldingCapacity) %>%
rename("tmin" = tmin_meanAnnAvg_CLIM,
"tmax" = tmax_meanAnnAvg_CLIM, #1
"tmean" = tmean_meanAnnAvg_CLIM,
"prcp" = prcp_meanAnnTotal_CLIM,
"t_warm" = T_warmestMonth_meanAnnAvg_CLIM,
"t_cold" = T_coldestMonth_meanAnnAvg_CLIM,
"prcp_wet" = precip_wettestMonth_meanAnnAvg_CLIM,
"prcp_dry" = precip_driestMonth_meanAnnAvg_CLIM,
"prcp_seasonality" = precip_Seasonality_meanAnnAvg_CLIM, #2
"prcpTempCorr" = PrecipTempCorr_meanAnnAvg_CLIM, #3
"abvFreezingMonth" = aboveFreezing_month_meanAnnAvg_CLIM,
"isothermality" = isothermality_meanAnnAvg_CLIM, #4
"annWatDef" = annWaterDeficit_meanAnnAvg_CLIM,
"annWetDegDays" = annWetDegDays_meanAnnAvg_CLIM,
"VPD_mean" = annVPD_mean_meanAnnAvg_CLIM,
"VPD_max" = annVPD_max_meanAnnAvg_CLIM, #5
"VPD_min" = annVPD_min_meanAnnAvg_CLIM, #6
"VPD_max_95" = annVPD_max_95percentile_CLIM,
"annWatDef_95" = annWaterDeficit_95percentile_CLIM,
"annWetDegDays_5" = annWetDegDays_5percentile_CLIM,
"frostFreeDays_5" = durationFrostFreeDays_5percentile_CLIM,
"frostFreeDays" = durationFrostFreeDays_meanAnnAvg_CLIM,
"soilDepth" = soilDepth, #7
"clay" = surfaceClay_perc,
"sand" = avgSandPerc_acrossDepth, #8
"coarse" = avgCoarsePerc_acrossDepth, #9
"carbon" = avgOrganicCarbonPerc_0_3cm, #10
"AWHC" = totalAvailableWaterHoldingCapacity,
## anomaly variables
tmean_anom = tmean_meanAnnAvg_3yrAnom, #15
tmin_anom = tmin_meanAnnAvg_3yrAnom, #16
tmax_anom = tmax_meanAnnAvg_3yrAnom, #17
prcp_anom = prcp_meanAnnTotal_3yrAnom, #18
t_warm_anom = T_warmestMonth_meanAnnAvg_3yrAnom, #19
t_cold_anom = T_coldestMonth_meanAnnAvg_3yrAnom, #20
prcp_wet_anom = precip_wettestMonth_meanAnnAvg_3yrAnom, #21
precp_dry_anom = precip_driestMonth_meanAnnAvg_3yrAnom, #22
prcp_seasonality_anom = precip_Seasonality_meanAnnAvg_3yrAnom, #23
prcpTempCorr_anom = PrecipTempCorr_meanAnnAvg_3yrAnom, #24
aboveFreezingMonth_anom = aboveFreezing_month_meanAnnAvg_3yrAnom, #25
isothermality_anom = isothermality_meanAnnAvg_3yrAnom, #26
annWatDef_anom = annWaterDeficit_meanAnnAvg_3yrAnom, #27
annWetDegDays_anom = annWetDegDays_meanAnnAvg_3yrAnom, #28
VPD_mean_anom = annVPD_mean_meanAnnAvg_3yrAnom, #29
VPD_min_anom = annVPD_min_meanAnnAvg_3yrAnom, #30
VPD_max_anom = annVPD_max_meanAnnAvg_3yrAnom, #31
VPD_max_95_anom = annVPD_max_95percentile_3yrAnom, #32
annWatDef_95_anom = annWaterDeficit_95percentile_3yrAnom, #33
annWetDegDays_5_anom = annWetDegDays_5percentile_3yrAnom , #34
frostFreeDays_5_anom = durationFrostFreeDays_5percentile_3yrAnom, #35
frostFreeDays_anom = durationFrostFreeDays_meanAnnAvg_3yrAnom #36
) %>%
dplyr::select(-c(tmin_meanAnnAvg_3yr:Start_3yr))
rm(climDat_temp)
gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 3083804 164.7 5090267 271.9 NA 3880161 207.3
## Vcells 950840264 7254.4 2664390240 20327.7 65536 2664106303 20325.6## now, scale the contemporary climate data for use in models
# get the scaling factors
scaleParams <- modDat_1_s %>%
#filter(Year == 2016) %>%
dplyr::select(tmin_s:AWHC_s) %>%
reframe(across(all_of(names(.)), attributes))
# apply the scaling factors to the contemporary climate data
namesToScale <- climDat %>%
dplyr::select(tmin:frostFreeDays, tmean_anom:frostFreeDays_anom, soilDepth:AWHC) %>%
names()
climDat_scaled <- map(namesToScale, .f = function(x) {
x_new <- (climDat[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
return(data.frame(x_new))
}) %>%
purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")
climDatPred <- climDat %>%
dplyr::select(NA_L1CODE:y) %>%
cbind(climDat_scaled)
names(climDatPred)[7:56] <- str_remove(names(climDatPred)[7:56], pattern = "_s$")
# apply the scaling factors to the contemporary climate data used for predictions over time
climDat_scaled <- map(namesToScale, .f = function(x) {
x_new <- (climDat_time[,x] - scaleParams[,paste0(x, "_s")]$`scaled:center`)/scaleParams[,paste0(x, "_s")]$`scaled:scale`
return(data.frame(x_new))
}) %>%
purrr::list_cbind()
names(climDat_scaled) <- paste0(namesToScale, "_s")
climDat_time_Pred <- climDat_time %>%
dplyr::select(NA_L1CODE:y,spatialID) %>%
cbind(climDat_scaled)
names(climDat_time_Pred)[9:58] <- str_remove(names(climDat_time_Pred)[9:58], pattern = "_s$")
rm(climDat_scaled)
gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 3091120 165.1 5090267 271.9 NA 5090267 271.9
## Vcells 976405202 7449.4 2664390240 20327.7 65536 2664106303 20325.6prednames_s <- modDat_1_s %>%
dplyr::select(tmin_s:AWHC_s) %>%
names()
prednames <- str_replace(prednames_s, pattern = "_s$", replacement = "")
## read in scaled data for future climate
forecastClimSoilsDatPred_1 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5.rds")
forecastClimSoilsDatPred_2 <- readRDS(file = "../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5.rds")
# read in unscaled data for future climate
forecastClimSoilsDat_1 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_BNU-ESM_rcp8_5_UNSCALED.rds")
forecastClimSoilsDat_2 <- readRDS("../../../Data_processed/CoverData/IntermediateAnalysisFiles/Final_ForecastedClimateDataAndSoilsDataForPredictions_IPSL-CM5A-MR_rcp8_5_UNSCALED.rds")First, we show predictions from the model of ecoregion classification
These maps show predictions of ecoregion classification using contemporary and modeled future climate data, where the predicted value is the probability that a location is forest. Note that while the climate data used in model predictions is either contemporary or from GCMs, the soils data used is static across time.
ecoMod <- readRDS("../../ecoregionClassification/ModelFitting/ModelObjects/EcoregionClassificationModel.rds")
## print out model statement
coefficientTable <- data.frame(coefficients(ecoMod))
responseVar <- "P(forest)"
coefficientTable$coefficientName <- rownames(coefficientTable)
coefficientTable$coefficients.ecoMod. <- round(coefficientTable$coefficients.ecoMod., 4)
# print out coefficients w/ coefficient names
tempNames <- paste0(
apply(coefficientTable, MARGIN = 1, FUN = function(x) {
if (x["coefficientName"] == "(Intercept)") {
paste0(x["coefficients.ecoMod."])
} else {
paste0(x["coefficients.ecoMod."], "*", x["coefficientName"])
}
}
),
collapse = " + "
)
# print the unscaled model statement
unscaledModelName <- paste0(responseVar, "~ 1/ (1 + exp(-(", tempNames,")))")
#print(unscaledModelName)# Here are maps of the ecoregion model predictions with both contemporary and modeled future climate data
## get model data used for ecoregion fitting
modDat_ecoregionFit <- readRDS("../../../Data_processed/EcoRegion_climSoilData.rds")
modDat_ecoregionFit <- modDat_ecoregionFit %>%
rename("Long" = x, "Lat" = y) %>%
dplyr::select(c(newRegion, #swe_meanAnnAvg_CLIM:
tmin_meanAnnAvg_CLIM:durationFrostFreeDays_meanAnnAvg_CLIM,
soilDepth , surfaceClay_perc ,
avgSandPerc_acrossDepth , avgCoarsePerc_acrossDepth,
avgOrganicCarbonPerc_0_3cm , totalAvailableWaterHoldingCapacity,
Long, Lat)) %>%
mutate(newRegion = as.factor(newRegion)) %>%
drop_na()
# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDatPred_unscaled <- climDat
names(climDatPred_unscaled)[c(5, 6, 7, 13, 10, 12, 98, 100, 101, 102)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
#rm(climDat)
#gc()
# predict with contemporary climate data
preds_byHand <- climDatPred_unscaled %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)
# predict with modeled future climate data v1
names(forecastClimSoilsDat_1)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future1_byHand <- forecastClimSoilsDat_1 %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
# predict with modeled future climate data v2
names(forecastClimSoilsDat_2)[c(8, 9, 10, 16, 13, 15, 48, 50, 51, 52)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
preds_future2_byHand <- forecastClimSoilsDat_2 %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
## Plot predictions with contemporary data
# make into a raster
plotObs <- preds_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
# get the extent of this particular raster, and crop it accordingly
tempExt <- crds(plotObs, na.rm = TRUE)
plotObs_2 <- plotObs %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_current <- ggplot() +
geom_spatraster(data = plotObs_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using
contemporary dayMet data")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)])
# for future v1
plotObs_future1 <- preds_future1_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
plotObs_future1_2 <- plotObs_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_future1 <- ggplot() +
geom_spatraster(data = plotObs_future1_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future1_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future
climate data from the BNU-ESM model")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_future1_2)[c(1,3)]) +
ylim(st_bbox(plotObs_future1_2)[c(2,4)])
# for future v2
plotObs_future2 <- preds_future2_byHand %>%
#drop_na(paste(response)) %>%
#slice_sample(n = 5e4) %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = "pred",
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
plotObs_future2_2 <- plotObs_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
map_ecoRegion_future2 <- ggplot() +
geom_spatraster(data = plotObs_future2_2) +
geom_sf(data = mapRegions, fill = NA, col = "orchid", lwd = .5) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>% st_crop(st_bbox(plotObs_future2_2)),fill=NA ) +
labs(title = paste0("Ecoregion classification predictions using future
climate data from the IPSL-CM5A-MR model")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey",
labels = c("grassShrub", "0.25", "0.50", "0.75", "forest")) +
xlim(st_bbox(plotObs_future2_2)[c(1,3)]) +
ylim(st_bbox(plotObs_future2_2)[c(2,4)])
## Plot the maps together
ggarrange(map_ecoRegion_current, map_ecoRegion_future1, map_ecoRegion_future2, nrow = 1) %>%
annotate_figure(fig.lab = "Model Predictions of Ecoregion Classification with Contemporary and Forecasted Climate and Data", fig.lab.size = 20,
bottom = text_grob("Note: Pink lines show the boundaries of contemporary grass/shrub and forest ecoregions, based on aggregated, EPA level 1 ecoregions"))
Now, we generate CONUS-wide model predictions of absolute cover
First, we use our models to predict absolute cover at the ecoregion or CONUS scale, depending on the functional type, using both contemporary climate and weather data, and then with future climate and weather data downscaled from the two GCMs. For some functional types we were able to fit a model that worked well across all of CONUS, while for other types, models performed better when we fit them to a single ecoregion at a time. Below is a list of the models we use:
- total herbaceous cover in the grass/shrub ecoregion
- total herbaceous cover in the forest ecoregion
- total tree cover in the grass/shrub ecoregion
- total tree cover in the forest ecoregion
- total shrub cover in CONUS
- total bare ground cover in CONUS
- The proportion of total tree cover that is needle-leaved in the grass/shrub ecoregion
- The proportion of total tree cover that is needle-leaved in the forest ecoregion
- The proportion of total tree cover that is broad-leaved in the grass/shrub ecoregion
- The proportion of total tree cover that is broad-leaved in the forest ecoregion
- The proportion of total herbaceous cover that is forbs in CONUS
- The proportion of total herbaceous cover that is C3 grass in CONUS
- The proportion of total herbaceous cover that is C4 grass in CONUS
*Note: Complete equations for each model used here can be found in the document “03_CompareModels.html”
# Total herbaceous cover:
# Grass/shrub: 1SE lambda model
totalHerb_GS <- readRDS("./models/betaLASSO/TotalHerbaceousCover_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Forest: best lambda model
totalHerb_F <- readRDS("./models/betaLASSO/TotalHerbaceousCover_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Total tree cover:
# Grass/shrub: best lambda model
totalTree_GS <- # include anomalies option #readRDS("./models/betaLASSO/TotalTreeCover_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
#remove anomalies option
readRDS("./models/betaLASSO/TotalTreeCover_shrubGrass_noTLP_FALSE_removeAnomaliesTRUE_bestLambdaGLM.rds")
# Forest: best lambda model
totalTree_F <- # include anomalies option
#remove anomalies option #
readRDS("./models/betaLASSO/TotalTreeCover_forest_noTLP_FALSE_removeAnomaliesTRUE_halfSELambdaGLM.rds")
# Bare ground:
# CONUS - 1SE lambda model
bareGround_CONUS <- readRDS("./models/betaLASSO/BareGroundCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# Shrub:
# CONUS - best lambda model
shrub_CONUS <- readRDS("./models/betaLASSO/ShrubCover_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Broad-leaved tree:
# Grass/shrub: best lambda model
BLtree_GS <- readRDS("./models/betaLASSO/AngioTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest: best lambda model
BLtree_F <- readRDS("./models/betaLASSO/AngioTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Needle-leaved tree:
# Grass/shrub: best lambda model
NLtree_GS <- readRDS("./models/betaLASSO/ConifTreeCover_prop_shrubGrass_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forest: best lambda model
NLtree_F <- readRDS("./models/betaLASSO/ConifTreeCover_prop_forest_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# Forb:
# CONUS: best lambda model
forb_CONUS <- readRDS("./models/betaLASSO/ForbCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
# C3 grass:
# CONUS: 1SE lambda model
C3grass_CONUS <- readRDS("./models/betaLASSO/C3GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_oneSELambdaGLM.rds")
# C4 grass:
# CONUS: best lambda model
C4grass_CONUS <- readRDS("./models/betaLASSO/C4GramCover_prop_CONUS_noTLP_FALSE_removeAnomaliesFALSE_trimAnom_bestLambdaGLM.rds")
## Now, predict absolute cover with scaled climate/weather/soils data
## with contemporary climate data
totalHerb_F_predsContemp <- makePredictions(predictionDF = climDatPred,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsContemp <- makePredictions(predictionDF = climDatPred,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsContemp<- makePredictions(predictionDF = climDatPred,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
contempPreds <- totalHerb_F_predsContemp %>%
cbind(totalHerb_GS_predsContemp %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsContemp %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsContemp %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsContemp %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsContemp %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsContemp %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsContemp %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsContemp %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsContemp %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsContemp %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsContemp %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsContemp %>% select(percForb_CONUS)) %>%
cbind(preds_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
# get predictions with first climate model
## with contemporary climate data
totalHerb_F_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture1 <- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsFuture1<- makePredictions(predictionDF = forecastClimSoilsDatPred_1,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsFuture1 <- totalHerb_F_predsFuture1 %>%
cbind(totalHerb_GS_predsFuture1 %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsFuture1 %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsFuture1 %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsFuture1 %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsFuture1 %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsFuture1 %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsFuture1 %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsFuture1 %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsFuture1 %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsFuture1 %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsFuture1 %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsFuture1 %>% select(percForb_CONUS)) %>%
cbind(preds_future1_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
# get predictions with second climate model
totalHerb_F_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsFuture2 <- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsFuture2<- makePredictions(predictionDF = forecastClimSoilsDatPred_2,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsFuture2 <- totalHerb_F_predsFuture2 %>%
cbind(totalHerb_GS_predsFuture2 %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsFuture2 %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsFuture2 %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsFuture2 %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsFuture2 %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsFuture2 %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsFuture2 %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsFuture2 %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsFuture2 %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsFuture2 %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsFuture2 %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsFuture2 %>% select(percForb_CONUS)) %>%
cbind(preds_future2_byHand %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)Once we generate predictions using each of these models, we transform the ecoregion-level predictions into CONUS-wide predictions by scaling them according to predicted ecoregion probability from the ecoregion classification model. e.g.:absolute total tree cover in CONUS = (absolute total tree cover in grass/shrub * P(grass/shrub)) + (absolute total tree cover in forest * P(forest))
## for contemporary data
contempPreds <- contempPreds %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)Then, for those models that are predicting proportions, we scale all of the relevant predictions so they sum to one. (e.g. Model predictions of the proportion of total herbaceous cover that is forbs, the proportion of total herbaceous cover that is C3 grass, and the proportion of total herbaceous cover that is C4 grass should sum to 1, so: scaled proportion of total herbaceous cover that is forbs = proportion of total herbaceous cover that is forbs/(proportion of total herbaceous cover that is forbs + proportion of total herbaceous cover that is C3 grass + proportion of total herbaceous cover that is C4 grass))
contempPreds <- contempPreds %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
predsFuture1 <- predsFuture1 %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
predsFuture2 <- predsFuture2 %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))Finally, we convert these predicted proportions into predictions of absolute cover. (e.g. absolute cover of needle-leaved trees = scaled proportion of total herbaceous cover that is forbs * absolute cover of total trees)
## for contemporary data
contempPreds <- contempPreds %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)Below are maps of predicted absolute cover for each functional group we modeled, after predictions have been scaled using the process outlined above. These maps show model-predicted cover using contemporary climate and weather data, residuals when comparing these contemporary predictions to observations, and changes in model predictions of cover when predictions are made using simulated future climate and weather data from two GCMs.
# prep observed data for use in figures
obsDat <- modDat_1_s %>%
select(Long, Lat, Year, ShrubCover, TotalHerbaceousCover, TotalTreeCover, BareGroundCover, C3GramCover_prop, C4GramCover_prop, ForbCover_prop, AngioTreeCover_prop, ConifTreeCover_prop) %>%
mutate(absNLTree_CONUS = ConifTreeCover_prop * TotalTreeCover,
absBLTree_CONUS = AngioTreeCover_prop * TotalTreeCover,
absC3grass_CONUS = C3GramCover_prop * TotalHerbaceousCover,
absC4grass_CONUS = C4GramCover_prop * TotalHerbaceousCover,
absForb_CONUS = ForbCover_prop * TotalHerbaceousCover
) %>%
rename(absShrub_CONUS = ShrubCover,
absTotalTree_CONUS = TotalTreeCover,
absTotalHerb_CONUS = TotalHerbaceousCover,
absBareGround_CONUS = BareGroundCover)
# list of absolute cover names
responseNames <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")
absoluteCoverMaps_contemp <- lapply(responseNames,
FUN = function(x) {
## contemporary absolute cover
# rasterize response
resp_rast <- contempPreds %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_2 <- resp_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# make a map
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of absolute cover of ",x,",
\n after scaling according to ecoregion prediction")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x)) +
theme_classic()
## calculate residuals
# make maps of residuals between absolute cover predictions and observations
# rasterize observations
plotObservations <- obsDat %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
plotObservations <- plotObservations %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate residuals (observed - predicted)
resids_rast <- plotObservations - resp_rast_2
## aggregate the residuals raster to make it easier to see
resids_rastCoarse <- resids_rast %>%
aggregate(fact = 2, fun = "mean", na.rm = TRUE)
# map the residuals
mapResids <- ggplot() +
geom_spatraster(data = resids_rastCoarse) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(resids_rast)),fill=NA ) +
labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-1,1)
) +
xlim(st_bbox(resids_rast)[c(1,3)]) +
ylim(st_bbox(resids_rast)[c(2,4)]) +
theme_classic()
# make a histogram of residuals
histResids <- ggplot() +
geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)),
col = "black", fill = "darkgrey") +
geom_vline(aes(xintercept = 0), linetype = 2) +
xlab(paste0("Residuals (obs - pred) ",x)) +
theme_classic()
## future model 1 absolute cover
# rasterize response
resp_rast_future1 <- predsFuture1 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future1_2 <- resp_rast_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
# make a map
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
## future model 2 absolute cover
# rasterize response
resp_rast_future2 <- predsFuture2 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future2_2 <- resp_rast_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
# make a map
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of absolute cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of absolute cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
return(list("rast" = resp_rast_2,
"map" = resp_map,
"hist" = resp_hist,
"rast_resids" = resids_rast,
"map_resids" = mapResids,
"hist_resids" = histResids,
"rast_future1" = resp_rast_future1_2,
"map_future1" = resp_map_future1,
"hist_future1" = resp_hist_future1,
"rast_future2" = resp_rast_future2_2,
"map_future2" = resp_map_future2,
"hist_future2" = resp_hist_future2))
})
names(absoluteCoverMaps_contemp) <- c("absTotalTree_CONUS", "absTotalHerb_CONUS", "absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS", "absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")ggarrange(
ggarrange(absoluteCoverMaps_contemp$absTotalTree_CONUS$map,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absTotalTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absTotalTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absTotalHerb_CONUS$map,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_resids,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future1,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$map_future2,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_resids,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future1,
absoluteCoverMaps_contemp$absTotalHerb_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absShrub_CONUS$map,
absoluteCoverMaps_contemp$absShrub_CONUS$map_resids,
absoluteCoverMaps_contemp$absShrub_CONUS$map_future1,
absoluteCoverMaps_contemp$absShrub_CONUS$map_future2,
absoluteCoverMaps_contemp$absShrub_CONUS$hist,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_resids,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_future1,
absoluteCoverMaps_contemp$absShrub_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absBareGround_CONUS$map,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_resids,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_future1,
absoluteCoverMaps_contemp$absBareGround_CONUS$map_future2,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_resids,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future1,
absoluteCoverMaps_contemp$absBareGround_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absNLTree_CONUS$map,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absNLTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absNLTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absBLTree_CONUS$map,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_resids,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_future1,
absoluteCoverMaps_contemp$absBLTree_CONUS$map_future2,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_resids,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future1,
absoluteCoverMaps_contemp$absBLTree_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absC3grass_CONUS$map,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_resids,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_future1,
absoluteCoverMaps_contemp$absC3grass_CONUS$map_future2,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_resids,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future1,
absoluteCoverMaps_contemp$absC3grass_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absC4grass_CONUS$map,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_resids,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_future1,
absoluteCoverMaps_contemp$absC4grass_CONUS$map_future2,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_resids,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future1,
absoluteCoverMaps_contemp$absC4grass_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(absoluteCoverMaps_contemp$absForb_CONUS$map,
absoluteCoverMaps_contemp$absForb_CONUS$map_resids,
absoluteCoverMaps_contemp$absForb_CONUS$map_future1,
absoluteCoverMaps_contemp$absForb_CONUS$map_future2,
absoluteCoverMaps_contemp$absForb_CONUS$hist,
absoluteCoverMaps_contemp$absForb_CONUS$hist_resids,
absoluteCoverMaps_contemp$absForb_CONUS$hist_future1,
absoluteCoverMaps_contemp$absForb_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ncol = 1
)
Now, we show quantile plots for model model predictions of absolute cover made with contemporary data. These plots compare final, scaled predictions of absolute cover to observed absolute cover. However, instead of showing raw data, which can be uninformative with large datasets such as ours, we show the average model predicted value or observed value within a given quantile of a given predictor variable. In this case, we are showing values within each percentile. We also show the inter-quantile range (IQR) of model-predicted and response values in each percentile, which appear as bands behind the points showing the mean values. Note that the predictor variables included in each quantile plot are the predictors that appear in all of the models used to generate a given absolute cover prediction (e.g. the predictors variables included in the quantile plot for absolute needle-leaved tree cover are the predictors that are included in any o fthe following models: total tree cover-grass/shrub, total tree cover-forest, prop. needle-leaved tree cover-grass/shrub, and prop. needle-leaved cover-forest). Note also that the predictor variables in these figures are scaled.
# predict ecoregion classification
modDat_ecoregionFitQuantile <- modDat_1_s %>%
#rename("Long" = x, "Lat" = y) %>%
dplyr::select(c(newRegion, t_warm, t_cold, prcp_wet, annWatDef, prcpTempCorr, isothermality, soilDepth, sand, coarse, carbon,
Long, Lat)) %>%
mutate(newRegion = as.factor(newRegion))
# get climate data from dayMet (d.f. is "climDatPred_unscaled")
modDat_ecoregionFitQuantile_unscaled <- modDat_ecoregionFitQuantile
names(modDat_ecoregionFitQuantile_unscaled)[c(2:11)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
# predict with contemporary climate data
preds_byHand_quantile <- modDat_ecoregionFitQuantile_unscaled %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
## now, predict for each functional type
# rename the data.frame so that the scaled predictors have the correct names for the models to predict with
scaledPrednames <- modDat_1_s %>%
select(tmin_s:AWHC_s) %>%
names()
modDat_quantile <- modDat_1_s %>%
select(ShrubCover:Lat, tmin_s:AWHC_s) %>%
rename_with(~str_remove(string = .x, pattern = "_s$"), any_of(scaledPrednames))
## with the observed data, convert proportions to absolute cover
modDat_quantile <- modDat_quantile %>%
mutate(C3GramCover_abs = C3GramCover_prop * TotalHerbaceousCover,
C4GramCover_abs = C4GramCover_prop * TotalHerbaceousCover,
ForbCover_abs = ForbCover_prop * TotalHerbaceousCover,
NLTreeCover_abs = ConifTreeCover_prop * TotalTreeCover,
BLTreeCover_abs = AngioTreeCover_prop * TotalTreeCover)
## with contemporary climate data
totalHerb_F_predsQuantile <- makePredictions(predictionDF = modDat_quantile,
modelObject = totalHerb_F, scale = FALSE) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsQuantile <- makePredictions(predictionDF = modDat_quantile,
modelObject = totalHerb_GS, scale = FALSE) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = totalTree_F, scale = FALSE) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = totalTree_GS, scale = FALSE) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = bareGround_CONUS, scale = FALSE) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = shrub_CONUS, scale = FALSE) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = BLtree_F, scale = FALSE) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = BLtree_GS, scale = FALSE) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = NLtree_F, scale = FALSE) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = NLtree_GS, scale = FALSE) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = C3grass_CONUS, scale = FALSE) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = C4grass_CONUS, scale = FALSE) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsQuantile<- makePredictions(predictionDF = modDat_quantile,
modelObject = forb_CONUS, scale = FALSE) %>%
rename(percForb_CONUS = "modelPreds")
preds_quantile <- totalHerb_F_predsQuantile %>%
cbind(totalHerb_GS_predsQuantile %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsQuantile %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsQuantile %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsQuantile %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsQuantile %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsQuantile %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsQuantile %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsQuantile %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsQuantile %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsQuantile %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsQuantile %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsQuantile %>% select(percForb_CONUS)) %>%
cbind(preds_byHand_quantile %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
# scale according to ecoregion
preds_quantile <- preds_quantile %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## For percentages (NL/BL tree and C4/C3/Forb), scale so that they sum to one
preds_quantile <- preds_quantile %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
## Convert percentages for level 2 cover into absolute cover values
preds_quantile <- preds_quantile %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)
## with the observed data, convert proportions to absolute cover
preds_quantile <- preds_quantile %>%
mutate(C3GramCover_abs = C3GramCover_prop * TotalHerbaceousCover,
C4GramCover_abs = C4GramCover_prop * TotalHerbaceousCover,
ForbCover_abs = ForbCover_prop * TotalHerbaceousCover,
NLTreeCover_abs = ConifTreeCover_prop * TotalTreeCover,
BLTreeCover_abs = AngioTreeCover_prop * TotalTreeCover)# total trees
# predictions for absolute needle leaved tree cover
deciles_totTree <- predvars2deciles(df = preds_quantile %>% select(TotalTreeCover, absTotalTree_CONUS, "tmean", "prcp", "prcp_dry", "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon" ,"coarse", "sand", "prcp_seasonality") %>% drop_na()
, response_vars = c("TotalTreeCover", "absTotalTree_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("tmean", "prcp", "prcp_dry", "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon" ,"coarse", "sand", "prcp_seasonality"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# NL trees
# predictions for absolute needle leaved tree cover
deciles_NLtree <- predvars2deciles(df = preds_quantile %>% select(NLTreeCover_abs, absNLTree_CONUS, tmean, prcp, prcp_dry , isothermality, AWHC, prcpTempCorr, clay, carbon, coarse, sand , prcp_seasonality, prcpTempCorr_anom , isothermality_anom) %>% drop_na()
, response_vars = c("NLTreeCover_abs", "absNLTree_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("tmean", "prcp", "prcp_dry" , "isothermality", "AWHC", "prcpTempCorr", "clay", "carbon", "coarse", "sand" , "prcp_seasonality", "prcpTempCorr_anom" , "isothermality_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc NL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# # BL trees
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_BLtree <- predvars2deciles(df = preds_quantile %>% select("BLTreeCover_abs", "absBLTree_CONUS", "tmean" ,"prcp" ,
"prcp_dry" , "isothermality" ,"AWHC" ,
"prcpTempCorr" , "clay", "carbon","coarse", "sand","prcp_seasonality",
"prcp_seasonality_anom","isothermality_anom", "prcpTempCorr_anom" ) %>% drop_na(),
response_vars = c("BLTreeCover_abs", "absBLTree_CONUS"), # name of observations, followed by name of predictions
pred_vars = c( "tmean" ,"prcp" ,
"prcp_dry" , "isothermality" ,"AWHC" ,
"prcpTempCorr" , "clay", "carbon","coarse", "sand","prcp_seasonality",
"prcp_seasonality_anom","isothermality_anom", "prcpTempCorr_anom" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# shrubs
# unique(c(names(coefficients(shrub_CONUS))))
deciles_shrub <- predvars2deciles(df = preds_quantile %>% select("ShrubCover", "absShrub_CONUS","prcp", "prcp_seasonality", "prcpTempCorr", "sand" ,"coarse", "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ) %>% drop_na(),
response_vars = c("ShrubCover", "absShrub_CONUS"), # name of observations, followed by name of predictions
pred_vars = c( "prcp", "prcp_seasonality", "prcpTempCorr", "sand" ,"coarse", "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ), # for predictors, combine list of predictors for all models
cut_points = seq(0, 1, 0.01)
)
# bare ground
# unique(c(names(coefficients())))
deciles_bareGround <- predvars2deciles(df = preds_quantile %>% select("BareGroundCover", "absBareGround_CONUS", "tmean" ,"prcpTempCorr" , "isothermality" , "annWetDegDays" ,
"sand" ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ) %>% drop_na(),
response_vars = c("BareGroundCover", "absBareGround_CONUS"), # name of observations, followed by name of predictions
pred_vars = c( "tmean" ,"prcpTempCorr" , "isothermality" , "annWetDegDays" ,
"sand" ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# total herbaceous cover
#unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS))))
deciles_totalHerb <- predvars2deciles(df = preds_quantile %>% select("TotalHerbaceousCover", "absTotalHerb_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean",
"prcp_anom", "clay", "carbon") %>% drop_na(),
response_vars = c("TotalHerbaceousCover", "absTotalHerb_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean",
"prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# forbs
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_forb <- predvars2deciles(df = preds_quantile %>% select("ForbCover_abs", "absForb_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay",
"carbon", "prcp_seasonality", "annWetDegDays" , "prcp_seasonality_anom", "prcpTempCorr_anom" , "annWetDegDays_anom") %>% drop_na(),
response_vars = c("ForbCover_abs", "absForb_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay",
"carbon", "prcp_seasonality", "annWetDegDays" , "prcp_seasonality_anom", "prcpTempCorr_anom" , "annWetDegDays_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# C3 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C3grass <- predvars2deciles(df = preds_quantile %>% select("C3GramCover_abs", "absC3grass_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC" , "isothermality_anom",
"tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom", "prcpTempCorr_anom", "annWetDegDays_anom", "prcp_seasonality" ) %>% drop_na(),
response_vars = c("C3GramCover_abs", "absC3grass_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC" , "isothermality_anom",
"tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom", "prcpTempCorr_anom", "annWetDegDays_anom", "prcp_seasonality" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# C4 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C4grass <- predvars2deciles(df = preds_quantile %>% select("C4GramCover_abs", "absC4grass_CONUS", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na(),
response_vars = c("C4GramCover_abs", "absC4grass_CONUS"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)# total trees trees
# absolute cover
quantPlot_totalTrees <- decile_dotplot_pq(df = deciles_totTree, response= c("TotalTreeCover", "absTotalTree_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Total Tree Cover: Quantile Plot")
quantPlot_totalTrees <- add_dotplot_inset(quantPlot_totalTrees, df = deciles_totTree, response = c("TotalTreeCover", "absTotalTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# needle-leaved trees
# absolute cover
quantPlot_NLtreeAbs <- decile_dotplot_pq(df = deciles_NLtree, response= c("NLTreeCover_abs", "absNLTree_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Needle-Leaved Tree Cover: Quantile Plot")
quantPlot_NLtreeAbs <- add_dotplot_inset(quantPlot_NLtreeAbs, df = deciles_NLtree, response = c("NLTreeCover_abs", "absNLTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# broad-leaved trees
# absolute cover
quantPlot_BLtreeAbs <- decile_dotplot_pq(df = deciles_BLtree, response= c("BLTreeCover_abs", "absBLTree_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Broad-Leaved Tree Cover: Quantile Plot")
quantPlot_BLtreeAbs <- add_dotplot_inset(quantPlot_BLtreeAbs, df = deciles_BLtree, response = c("BLTreeCover_abs", "absBLTree_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# shrubs
quantPlot_shrubs <- decile_dotplot_pq(df = deciles_shrub, response= c("ShrubCover", "absShrub_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Shrub Cover: Quantile Plot")
quantPlot_shrubs <- add_dotplot_inset(quantPlot_shrubs, df = deciles_shrub, response = c("ShrubCover", "absShrub_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# bare ground
quantPlot_bareGround <- decile_dotplot_pq(df = deciles_bareGround, response= c("BareGroundCover", "absBareGround_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Bare Ground Cover: Quantile Plot")
quantPlot_bareGround <- add_dotplot_inset(quantPlot_bareGround, df = deciles_bareGround, response = c("BareGroundCover", "absBareGround_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# total herbaceous
quantPlot_totHerb <- decile_dotplot_pq(df = deciles_totalHerb, response= c("TotalHerbaceousCover", "absTotalHerb_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Total Herbaceous Cover: Quantile Plot")
quantPlot_totHerb <- add_dotplot_inset(quantPlot_totHerb, df = deciles_totalHerb, response = c("TotalHerbaceousCover", "absTotalHerb_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# forbs
quantPlot_forbs <- decile_dotplot_pq(df = deciles_forb, response= c("ForbCover_abs", "absForb_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Forb Cover: Quantile Plot")
quantPlot_forbs <- add_dotplot_inset(quantPlot_forbs, df = deciles_forb, response = c("ForbCover_abs", "absForb_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# C3 grass
quantPlot_c3grass <- decile_dotplot_pq(df = deciles_C3grass, response= c("C3GramCover_abs", "absC3grass_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute C3 Grass Cover: Quantile Plot")
quantPlot_c3grass <- add_dotplot_inset(quantPlot_c3grass, df = deciles_C3grass, response = c("C3GramCover_abs", "absC3grass_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
#C4 grass
quantPlot_c4grass <- decile_dotplot_pq(df = deciles_C4grass, response= c("C4GramCover_abs", "absC4grass_CONUS"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute C4 Grass Cover: Quantile Plot")
quantPlot_c4grass <- add_dotplot_inset(quantPlot_c4grass, df = deciles_C4grass, response = c("C4GramCover_abs", "absC4grass_CONUS"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)## now, display the figures
quantPlot_totalTrees
quantPlot_totHerb
quantPlot_shrubs
quantPlot_bareGround
quantPlot_NLtreeAbs
quantPlot_BLtreeAbs
quantPlot_c3grass
quantPlot_c4grass
quantPlot_forbs
Below, we show the RMSE and bias of the final, scaled estimates of absolute cover for each functional type in comparison to the observed data.
long_absPreds <- preds_quantile %>%
mutate(uniqueID = 1:nrow(preds_quantile)) %>%
select(uniqueID, absTotalHerb_CONUS, absTotalTree_CONUS, absShrub_CONUS, absBareGround_CONUS,
absNLTree_CONUS, absBLTree_CONUS, absC3grass_CONUS, absC4grass_CONUS, absForb_CONUS) %>%
rename(TotalHerb = absTotalHerb_CONUS,
TotalTree = absTotalTree_CONUS,
Shrub = absShrub_CONUS,
BareGround = absBareGround_CONUS,
NLTree = absNLTree_CONUS,
BLTree = absBLTree_CONUS,
C3gram = absC3grass_CONUS,
C4gram = absC4grass_CONUS,
Forb = absForb_CONUS) %>%
pivot_longer(cols = c(TotalHerb, TotalTree, Shrub, BareGround, NLTree, BLTree, C3gram, C4gram, Forb),
names_to = "names",
values_to = "preds_values"
) %>%
select(uniqueID, names, preds_values)
long_absObs <- preds_quantile %>%
mutate(uniqueID = 1:nrow(preds_quantile)) %>%
select(uniqueID, TotalHerbaceousCover, TotalTreeCover, ShrubCover, BareGroundCover,
NLTreeCover_abs, BLTreeCover_abs, C3GramCover_abs, C4GramCover_abs, ForbCover_abs) %>%
rename(TotalHerb = TotalHerbaceousCover,
TotalTree = TotalTreeCover,
Shrub = ShrubCover,
BareGround = BareGroundCover,
NLTree = NLTreeCover_abs,
BLTree = BLTreeCover_abs,
C3gram = C3GramCover_abs,
C4gram = C4GramCover_abs,
Forb = ForbCover_abs) %>%
pivot_longer(cols = c(TotalHerb, TotalTree, Shrub, BareGround, NLTree, BLTree, C3gram, C4gram, Forb),
names_to = "names",
values_to = "obs_values"
) %>%
select(uniqueID, names, obs_values)
longObsPreds <- long_absObs %>%
left_join(long_absPreds)
## calculate bias
(summaryTableAbs <- longObsPreds %>%
mutate(resids = obs_values - preds_values) %>%
group_by(names) %>%
dplyr::summarise(bias = mean(resids, na.rm = TRUE),
RMSE = round(sqrt(mean(resids^2, na.rm = TRUE)),4)
) %>%
slice(8, 9, 7, 2, 6, 1, 3, 4, 5) %>% kable(
format = "html",
col.names = c("Model Name", "Model bias mean(obs. - pred.)", "model RMSE")
))| Model Name | Model bias mean(obs. - pred.) | model RMSE |
|---|---|---|
| TotalHerb | 0.0020453 | 0.2181 |
| TotalTree | -0.0064183 | 0.1815 |
| Shrub | 0.0000475 | 0.1423 |
| BareGround | 0.0000404 | 0.1327 |
| NLTree | 0.0698285 | 0.2010 |
| BLTree | 0.0307924 | 0.1589 |
| C3gram | 0.0319504 | 0.1529 |
| C4gram | 0.0059215 | 0.0846 |
| Forb | -0.0041207 | 0.1314 |
Finally, we show how these model predictions made using contemporary climate and weather data change over time from 2010 to 2020.
# get climate data from dayMet (d.f. is "climDatPred_unscaled")
climDat_time_unscaled <- climDat_time
names(climDat_time_unscaled)[c(5, 6, 7, 13, 10, 12, 99, 101, 102, 103)] <- c("T_warmestMonth_meanAnnAvg_CLIM", "T_coldestMonth_meanAnnAvg_CLIM", "precip_wettestMonth_meanAnnAvg_CLIM", "annWaterDeficit_meanAnnAvg_CLIM", "PrecipTempCorr_meanAnnAvg_CLIM", "isothermality_meanAnnAvg_CLIM", "soilDepth", "avgSandPerc_acrossDepth", "avgCoarsePerc_acrossDepth", "avgOrganicCarbonPerc_0_3cm")
# predict with contemporary climate data
preds_byHand_time <- climDat_time_unscaled %>%
mutate(pred = 1/(1 + exp(-( 9.8726 + -0.2999*T_warmestMonth_meanAnnAvg_CLIM + 0.2456*T_coldestMonth_meanAnnAvg_CLIM + 0.0106*precip_wettestMonth_meanAnnAvg_CLIM + -0.0621*annWaterDeficit_meanAnnAvg_CLIM + -2.7863*PrecipTempCorr_meanAnnAvg_CLIM + 0.0540*isothermality_meanAnnAvg_CLIM + -0.0076*soilDepth + 0.0335*avgSandPerc_acrossDepth + 0.0310*avgCoarsePerc_acrossDepth + 0.2726*avgOrganicCarbonPerc_0_3cm))))
#predictionsModel <- predict(object = ecoMod, type = "response", newdata = climDatPred_unscaled)
## with contemporary climate data
totalHerb_F_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalHerb_F) %>%
rename(absTotalHerb_F = "modelPreds")
totalHerb_GS_predsOverTime <- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalHerb_GS) %>%
rename(absTotalHerb_GS = "modelPreds")
totalTree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalTree_F) %>%
rename(absTotalTree_F = "modelPreds")
totalTree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = totalTree_GS) %>%
rename(absTotalTree_GS = "modelPreds")
bareGround_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = bareGround_CONUS) %>%
rename(absBareGround_CONUS = "modelPreds")
shrub_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = shrub_CONUS) %>%
rename(absShrub_CONUS = "modelPreds")
BLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = BLtree_F) %>%
rename(percBLTree_F = "modelPreds")
BLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = BLtree_GS) %>%
rename(percBLTree_GS = "modelPreds")
NLtree_F_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = NLtree_F) %>%
rename(percNLTree_F = "modelPreds")
NLtree_GS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = NLtree_GS) %>%
rename(percNLTree_GS = "modelPreds")
C3grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = C3grass_CONUS) %>%
rename(percC3grass_CONUS = "modelPreds")
C4grass_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = C4grass_CONUS) %>%
rename(percC4grass_CONUS = "modelPreds")
forb_CONUS_predsOverTime<- makePredictions(predictionDF = climDat_time_Pred,
modelObject = forb_CONUS) %>%
rename(percForb_CONUS = "modelPreds")
predsOverTime <- totalHerb_F_predsOverTime %>%
cbind(totalHerb_GS_predsOverTime %>% select(absTotalHerb_GS)) %>%
cbind(totalTree_F_predsOverTime %>% select(absTotalTree_F)) %>%
cbind(totalTree_GS_predsOverTime %>% select(absTotalTree_GS)) %>%
cbind(bareGround_CONUS_predsOverTime %>% select(absBareGround_CONUS)) %>%
cbind(shrub_CONUS_predsOverTime %>% select(absShrub_CONUS)) %>%
cbind(BLtree_F_predsOverTime %>% select(percBLTree_F)) %>%
cbind(BLtree_GS_predsOverTime %>% select(percBLTree_GS)) %>%
cbind(NLtree_F_predsOverTime %>% select(percNLTree_F)) %>%
cbind(NLtree_GS_predsOverTime %>% select(percNLTree_GS)) %>%
cbind(C3grass_CONUS_predsOverTime %>% select(percC3grass_CONUS)) %>%
cbind(C4grass_CONUS_predsOverTime %>% select(percC4grass_CONUS)) %>%
cbind(forb_CONUS_predsOverTime %>% select(percForb_CONUS)) %>%
cbind(preds_byHand_time %>% select(pred)) %>%
rename(prob_Forest = "pred") %>%
mutate(prob_grassShrub = 1 - prob_Forest)
## now, scale predictions according to ecoregion percentage
predsOverTime <- predsOverTime %>%
mutate(absTotalHerb_CONUS = (prob_grassShrub * absTotalHerb_GS) + (prob_Forest * absTotalHerb_F),
absTotalTree_CONUS = (prob_grassShrub * absTotalTree_GS) + (prob_Forest * absTotalTree_F),
percBLTree_CONUS = (prob_grassShrub * percBLTree_GS) + (prob_Forest * percBLTree_F),
percNLTree_CONUS = (prob_grassShrub * percNLTree_GS) + (prob_Forest * percNLTree_F)
)
## then, scale the percentages so they sum to 1
predsOverTime <- predsOverTime %>%
mutate(percBLTree_scaled_CONUS = percBLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percNLTree_scaled_CONUS = percNLTree_CONUS/(percBLTree_CONUS + percNLTree_CONUS),
percC3grass_scaled_CONUS = percC3grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percC4grass_scaled_CONUS = percC4grass_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS),
percForb_scaled_CONUS = percForb_CONUS/(percC3grass_CONUS + percC4grass_CONUS + percForb_CONUS))
## convert percentages into absolute cover for level 2 cover classes
predsOverTime <- predsOverTime %>%
mutate(absNLTree_CONUS = (percNLTree_scaled_CONUS * absTotalTree_CONUS),
absBLTree_CONUS = (percBLTree_scaled_CONUS * absTotalTree_CONUS),
absC3grass_CONUS = (percC3grass_scaled_CONUS * absTotalHerb_CONUS),
absC4grass_CONUS = (percC4grass_scaled_CONUS * absTotalHerb_CONUS),
absForb_CONUS = (percForb_scaled_CONUS * absTotalHerb_CONUS)
)# change spatialIDs to discrete numbers, since the long form is causing some sort of floating point issue?
spatID_lu <- data.frame("spatialID" = unique(predsOverTime$spatialID),
"uniqueID" = 1:length(predsOverTime$spatialID))
predsOverTime <- predsOverTime %>%
left_join(spatID_lu)
# make into a long format
predsOverTime_long <- predsOverTime %>%
select(year,x, y, spatialID, uniqueID, absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS) %>%
pivot_longer(cols = c(absShrub_CONUS, absBareGround_CONUS, absTotalHerb_CONUS, absTotalTree_CONUS, absNLTree_CONUS:absForb_CONUS)) %>%
mutate(year = as.integer(year))
# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>%
group_by(uniqueID, name) %>%
summarize(cv = sd(value)/mean(value))
## add x,y
predsOverTime_cv <- predsOverTime_cv %>%
left_join(predsOverTime %>% select(uniqueID, x, y))
# visualize change over time for a selection of points
coverOverTimeForSelectLocations <- predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("absShrub_CONUS", "absBareGround_CONUS", "absNLTree_CONUS",
"absBLTree_CONUS", "absC3grass_CONUS", "absC4grass_CONUS", "absForb_CONUS")) %>%
ggplot() +
facet_wrap(~uniqueID) +
geom_area(aes(x = year, y = value, fill = name)) +
geom_line(data = predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("absTotalHerb_CONUS", "absTotalTree_CONUS")),
aes(x = year, y = value, linetype = name), col = "black") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = "Change in modeled absolute cover by functional type from \n 2010 to 2020 for a subset of locations")
mapOfSelectLocations <- predsOverTime_long %>%
select(x, y, uniqueID) %>%
filter(uniqueID %in% 1:50) %>%
unique() %>%
ggplot() +
geom_sf(data = cropped_states, aes()) +
geom_point(aes(x, y)) +
geom_label_repel(aes(x, y, label = uniqueID)) +
labs(title = "Locations in above plot")
ggarrange(coverOverTimeForSelectLocations, mapOfSelectLocations, ncol = 1, heights = c(1,.5))
predsOverTime_cv %>%
ggplot() +
geom_boxplot(aes(x = name, y = cv, col = name)) +
theme_minimal() +
labs(title = "Coefficients of variation of model-predicted absolute \n cover at 1000 random locations for each functional group from 2010-2020",
x = "functional group",
y = "coefficient of variation") +
theme(axis.text.x = element_text(angle = 90))
# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>%
# rename(geometry = Shape)
ggplot() +
facet_wrap(~name) +
geom_sf(data = cropped_states, aes()) +
stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) +
guides(fill = guide_legend(title = "CV")) +
theme_minimal() +
labs(title = "The spatial distribution of the coefficients of variation of \nmodel-predicted cover from 2010-2020 for 1000 random locations",
subtitle = "Note that values are aggregated accross space to enhance readability")
Now, we generate CONUS-wide model predictions of relative cover
We use the model predictions of absolute cover we generated above, which are CONUS-wide and have been scaled according to predicted ecoregions, to calculate relative cover for each functional group. We consider the overstory (trees) and understory (all other functional types) separately in our relativization scheme, such that the values of relative cover are identical to the values of absolute cover for total trees, needle-leaved trees, and broad-leaved trees. For the understory (comprised of all herbaceous functional types, shrubs, and bare ground), total relative cover must sum to 1, so absolute cover for these functional groups is scaled accordingly.
## for contemporary data
contempPreds <- contempPreds %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
## for first climate model
predsFuture1 <- predsFuture1 %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
## for second climate model
predsFuture2 <- predsFuture2 %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)Below are maps of predicted relative cover for each functional group we modeled. These maps show model-predicted relative cover using contemporary climate and weather data, residuals when comparing these contemporary predictions to observations, and changes in model predictions of cover when predictions are made using simulated future climate and weather data from two GCMs. This figure also includes a map of the total understory absolute cover, which was the scaling factor used to calculate relative cover for each understory group. Note that relative cover maps for total trees, needle-leaved trees, and broad-leaved trees are not shown, since, for those functional groups, relative cover is identical to absolute cover.
# list of absolute cover names
responseNames <- c(#"relCoverB_totalTree",
"relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
#"relCoverB_NLTree", "relCoverB_BLTree",
"relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
"total_NonTree_AbsCover_CONUS")
# prepare the observed data for raster calculations (calculate relative observed cover)
obsDat <- obsDat %>%
mutate(
total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + absShrub_CONUS + absBareGround_CONUS) # for observations
) %>%
mutate(
# making observations relative
relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS,
)
relCover_B_Maps_contemp <- lapply(responseNames,
FUN = function(x) {
## contemporary absolute cover
# rasterize response
resp_rast <- contempPreds %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_2 <- resp_rast %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
## now, rasterize the observed data and calculate residualks
plotObservations <- obsDat %>%
terra::vect(geom = c("Long", "Lat")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
plotObservations <- plotObservations %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate residuals (observed - predicted)
resids_rast <- plotObservations - resp_rast_2
## aggregate the residuals raster to make it easier to see
resids_rastCoarse <- resids_rast %>%
aggregate(fact = 2, fun = "mean", na.rm = TRUE)
# map the residuals
mapResids <- ggplot() +
geom_spatraster(data = resids_rastCoarse) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(resids_rast)),fill=NA ) +
labs(title = paste0("Residuals (observations - predicted absolute cover) for ",x)) +
scale_fill_gradient2(low = "red",
mid = "white" ,
high = "blue" ,
midpoint = 0, na.value = "lightgrey",
limits = c(-1,1)
) +
xlim(st_bbox(resids_rast)[c(1,3)]) +
ylim(st_bbox(resids_rast)[c(2,4)]) +
theme_classic()
# make a histogram of residuals
histResids <- ggplot() +
geom_density(aes(terra::values(resids_rast$mean, na.rm = TRUE, mat = FALSE)),
col = "black", fill = "darkgrey") +
geom_vline(aes(xintercept = 0), linetype = 2) +
xlab(paste0("Residuals (obs - pred) ",x)) +
theme_classic()
## future model 1 absolute cover
# rasterize response
resp_rast_future1 <- predsFuture1 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future1_2 <- resp_rast_future1 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future1 <- resp_rast_future1_2 - resp_rast_2
## future model 2 absolute cover
# rasterize response
resp_rast_future2 <- predsFuture2 %>%
terra::vect(geom = c("x", "y")) %>%
terra::set.crs(crs(test_rast)) %>%
terra::rasterize(y = test_rast,
field = x,
fun = mean, na.rm = TRUE)
# get the extent of this particular raster, and crop it accordingly
resp_rast_future2_2 <- resp_rast_future2 %>%
crop(ext(min(tempExt[,1]), max(tempExt[,1]),
min(tempExt[,2]), max(tempExt[,2]))
)
# calculate differences between contemporary preds and future 1 preds (deltas)
rast_deltas_future2 <- resp_rast_future2_2 - resp_rast_2
if (x == "total_NonTree_AbsCover_CONUS") {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
scale_fill_gradient2(low = "lightblue",
#mid = "wheat" ,
high = "darkblue" ,
#midpoint = 0,
limits = c(0,max(terra::values(resp_rast_2$mean), na.rm= TRUE)),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of Total Absolute Cover (non-tree) with contemporary data")) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Total Absolute Cover (non-tree) \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of Cover (non-tree), \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-2,2), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
} else {
# make a map of predictions w/ current data
resp_map <- ggplot() +
geom_spatraster(data = resp_rast_2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Predictions of relative cover of ",x,
", \n when all functional groups except trees sum to 1")) +
scale_fill_gradient2(low = "brown",
mid = "wheat" ,
high = "darkgreen" ,
midpoint = 0, limits = c(0,1),
na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist <- ggplot() +
geom_density(aes(contempPreds[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x)) +
theme_classic()
# make a map of future preds w/ model 1
resp_map_future1 <- ggplot() +
geom_spatraster(data = rast_deltas_future1) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",x,", \n compared to contemporary predictions (future - contemprary) \n - climate model: BNU-ESM")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future1 <- ggplot() +
geom_density(aes(predsFuture1[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = BNU-ESM)")) +
theme_classic()
# make a map of future preds w/ model 2
resp_map_future2 <- ggplot() +
geom_spatraster(data = rast_deltas_future2) +
geom_sf(data=cropped_states %>% st_transform(crs = st_crs(test_rast)) %>%
st_crop(st_bbox(plotObs_2)),fill=NA ) +
labs(title = paste0("Future predictions of relative cover of ",
x,
", \n compared to contemporary predictions (future - contemprary) \n- climate model: IPSL-CM5A-MR (France)")) +
scale_fill_gradient2(low = "orange",
mid = "white" ,
high = "purple" ,
midpoint = 0, limits = c(-1,1), na.value = "lightgrey") +
xlim(st_bbox(plotObs_2)[c(1,3)]) +
ylim(st_bbox(plotObs_2)[c(2,4)]) +
theme_classic()
# make a histogram
resp_hist_future2 <- ggplot() +
geom_density(aes(predsFuture2[,x]), col = "black", fill = "darkgrey") +
xlab(paste0("Predictions of relative cover of ",x,
" \n in future scenario (model = IPSL-CM5A-MR (France))")) +
theme_classic()
}
return(list("rast" = resp_rast_2,
"map" = resp_map,
"hist" = resp_hist,
"rast_resids" = resids_rast,
"map_resids" = mapResids,
"hist_resids" = histResids,
"rast_future1" = resp_rast_future1_2,
"map_future1" = resp_map_future1,
"hist_future1" = resp_hist_future1,
"rast_future2" = resp_rast_future2_2,
"map_future2" = resp_map_future2,
"hist_future2" = resp_hist_future2))
})
names(relCover_B_Maps_contemp) <- c(#"relCoverB_totalTree",
"relCoverB_totalHerb", "relCoverB_shrub", "relCoverB_bareGround",
#"relCoverB_NLTree", "relCoverB_BLTree",
"relCoverB_C3Grass", "relCoverB_C4Grass", "relCoverB_Forb",
"total_NonTree_AbsCover_CONUS")ggarrange(
# ggarrange(relCover_B_Maps_contemp$relCoverB_totalTree$map,
# relCover_B_Maps_contemp$relCoverB_totalTree$map_resids,
# relCover_B_Maps_contemp$relCoverB_totalTree$map_future1,
# relCover_B_Maps_contemp$relCoverB_totalTree$map_future2,
# relCover_B_Maps_contemp$relCoverB_totalTree$hist,
# relCover_B_Maps_contemp$relCoverB_totalTree$hist_resids,
# relCover_B_Maps_contemp$relCoverB_totalTree$hist_future1,
# relCover_B_Maps_contemp$relCoverB_totalTree$hist_future2,
# ncol = 4,
# nrow = 2,
# heights = c(1,.35)
# ),
ggarrange(relCover_B_Maps_contemp$relCoverB_totalHerb$map,
relCover_B_Maps_contemp$relCoverB_totalHerb$map_resids,
relCover_B_Maps_contemp$relCoverB_totalHerb$map_future1,
relCover_B_Maps_contemp$relCoverB_totalHerb$map_future2,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist_resids,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future1,
relCover_B_Maps_contemp$relCoverB_totalHerb$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_shrub$map,
relCover_B_Maps_contemp$relCoverB_shrub$map_resids,
relCover_B_Maps_contemp$relCoverB_shrub$map_future1,
relCover_B_Maps_contemp$relCoverB_shrub$map_future2,
relCover_B_Maps_contemp$relCoverB_shrub$hist,
relCover_B_Maps_contemp$relCoverB_shrub$hist_resids,
relCover_B_Maps_contemp$relCoverB_shrub$hist_future1,
relCover_B_Maps_contemp$relCoverB_shrub$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_bareGround$map,
relCover_B_Maps_contemp$relCoverB_bareGround$map_resids,
relCover_B_Maps_contemp$relCoverB_bareGround$map_future1,
relCover_B_Maps_contemp$relCoverB_bareGround$map_future2,
relCover_B_Maps_contemp$relCoverB_bareGround$hist,
relCover_B_Maps_contemp$relCoverB_bareGround$hist_resids,
relCover_B_Maps_contemp$relCoverB_bareGround$hist_future1,
relCover_B_Maps_contemp$relCoverB_bareGround$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
# ggarrange(relCover_B_Maps_contemp$relCoverB_NLTree$map,
# relCover_B_Maps_contemp$relCoverB_NLTree$map_resids,
# relCover_B_Maps_contemp$relCoverB_NLTree$map_future1,
# relCover_B_Maps_contemp$relCoverB_NLTree$map_future2,
# relCover_B_Maps_contemp$relCoverB_NLTree$hist,
# relCover_B_Maps_contemp$relCoverB_NLTree$hist_resids,
# relCover_B_Maps_contemp$relCoverB_NLTree$hist_future1,
# relCover_B_Maps_contemp$relCoverB_NLTree$hist_future2,
# ncol = 4,
# nrow = 2,
# heights = c(1,.35)
# ),
# ggarrange(relCover_B_Maps_contemp$relCoverB_BLTree$map,
# relCover_B_Maps_contemp$relCoverB_BLTree$map_resids,
# relCover_B_Maps_contemp$relCoverB_BLTree$map_future1,
# relCover_B_Maps_contemp$relCoverB_BLTree$map_future2,
# relCover_B_Maps_contemp$relCoverB_BLTree$hist,
# relCover_B_Maps_contemp$relCoverB_BLTree$hist_resids,
# relCover_B_Maps_contemp$relCoverB_BLTree$hist_future1,
# relCover_B_Maps_contemp$relCoverB_BLTree$hist_future2,
# ncol = 4,
# nrow = 2,
# heights = c(1,.35)
# ),
ggarrange(relCover_B_Maps_contemp$relCoverB_C3Grass$map,
relCover_B_Maps_contemp$relCoverB_C3Grass$map_resids,
relCover_B_Maps_contemp$relCoverB_C3Grass$map_future1,
relCover_B_Maps_contemp$relCoverB_C3Grass$map_future2,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist_resids,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future1,
relCover_B_Maps_contemp$relCoverB_C3Grass$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_C4Grass$map,
relCover_B_Maps_contemp$relCoverB_C4Grass$map_resids,
relCover_B_Maps_contemp$relCoverB_C4Grass$map_future1,
relCover_B_Maps_contemp$relCoverB_C4Grass$map_future2,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist_resids,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future1,
relCover_B_Maps_contemp$relCoverB_C4Grass$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$relCoverB_Forb$map,
relCover_B_Maps_contemp$relCoverB_Forb$map_resids,
relCover_B_Maps_contemp$relCoverB_Forb$map_future1,
relCover_B_Maps_contemp$relCoverB_Forb$map_future2,
relCover_B_Maps_contemp$relCoverB_Forb$hist,
relCover_B_Maps_contemp$relCoverB_Forb$hist_resids,
relCover_B_Maps_contemp$relCoverB_Forb$hist_future1,
relCover_B_Maps_contemp$relCoverB_Forb$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ggarrange(relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_resids,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future1,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$map_future2,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_resids,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future1,
relCover_B_Maps_contemp$total_NonTree_AbsCover_CONUS$hist_future2,
ncol = 4,
nrow = 2,
heights = c(1,.35)
),
ncol = 1
)
Now, we show quantile plots for model predictions of relative cover made with contemporary data. These plots compare final predictions of relative cover to observed relative cover. Quantile plots for relative cover of tree functional types are now show, since they are identical to those for absolute cover.
# First, we need to relativize the predictions and the observations (stored in preds_quantiles)
preds_quantile <- preds_quantile %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS), # for model predictions
total_NonTree_AbsCover_obs = (TotalHerbaceousCover + ShrubCover + BareGroundCover) # for observations
) %>%
# making model predictions relative
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS,
# making observations relative
relCoverB_totalTree_Obs = TotalTreeCover,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb_Obs = TotalHerbaceousCover/total_NonTree_AbsCover_obs,
relCoverB_shrub_Obs = ShrubCover/total_NonTree_AbsCover_obs,
relCoverB_bareGround_Obs = BareGroundCover/total_NonTree_AbsCover_obs,
relCoverB_NLTree_Obs = NLTreeCover_abs,
relCoverB_BLTree_Obs = BLTreeCover_abs,
relCoverB_C3Grass_Obs = C3GramCover_abs/total_NonTree_AbsCover_obs,
relCoverB_C4Grass_Obs = C4GramCover_abs/total_NonTree_AbsCover_obs,
relCoverB_Forb_Obs = ForbCover_abs/total_NonTree_AbsCover_obs,
)## exclude trees, since the predictions are the same as the absolute cover shown above
# shrubs
# unique(c(names(coefficients(shrub_CONUS))))
deciles_shrub_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_shrub_Obs", "relCoverB_shrub","prcp", "prcp_seasonality", "prcpTempCorr", "sand" ,"coarse", "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ) %>% drop_na(),
response_vars = c("relCoverB_shrub_Obs", "relCoverB_shrub"), # name of observations, followed by name of predictions
pred_vars = c( "prcp", "prcp_seasonality", "prcpTempCorr", "sand" ,"coarse", "isothermality_anom", "AWHC" ,"isothermality", "annWetDegDays", "prcp_seasonality_anom" ,"tmean" ), # for predictors, combine list of predictors for all models
cut_points = seq(0, 1, 0.01)
)
# bare ground
# unique(c(names(coefficients())))
deciles_bareGround_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_bareGround_Obs", "relCoverB_bareGround", "tmean" ,"prcpTempCorr" , "isothermality" , "annWetDegDays" ,
"sand" ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ) %>% drop_na(),
response_vars = c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), # name of observations, followed by name of predictions
pred_vars = c( "tmean" ,"prcpTempCorr" , "isothermality" , "annWetDegDays" ,
"sand" ,"coarse" , "AWHC", "isothermality_anom" ,"annWetDegDays_anom" ,"prcp" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# total herbaceous cover
#unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS))))
deciles_totalHerb_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_totalHerb_Obs", "relCoverB_totalHerb", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean",
"prcp_anom", "clay", "carbon") %>% drop_na(),
response_vars = c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean",
"prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# forbs
# unique(c(names(coefficients(totalTree_F)), names(coefficients(totalTree_GS)), names(coefficients(BLtree_GS)), names(coefficients(BLtree_F))))
deciles_forb_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_Forb_Obs", "relCoverB_Forb", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay",
"carbon", "prcp_seasonality", "annWetDegDays" , "prcp_seasonality_anom", "prcpTempCorr_anom" , "annWetDegDays_anom") %>% drop_na(),
response_vars = c("relCoverB_Forb_Obs", "relCoverB_Forb"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean" , "prcp_anom", "clay",
"carbon", "prcp_seasonality", "annWetDegDays" , "prcp_seasonality_anom", "prcpTempCorr_anom" , "annWetDegDays_anom"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# C3 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C3grass_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_C3Grass_Obs", "relCoverB_C3Grass", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC" , "isothermality_anom",
"tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom", "prcpTempCorr_anom", "annWetDegDays_anom", "prcp_seasonality" ) %>% drop_na(),
response_vars = c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC" , "isothermality_anom",
"tmean", "prcp_anom", "clay", "carbon", "annWetDegDays", "prcp_seasonality_anom", "prcpTempCorr_anom", "annWetDegDays_anom", "prcp_seasonality" ), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)
# C4 grass
# unique(c(names(coefficients(totalHerb_F)), names(coefficients(totalHerb_GS)), names(coefficients(C3grass_CONUS))))
deciles_C4grass_Obs <- predvars2deciles(df = preds_quantile %>% select("relCoverB_C4Grass_Obs", "relCoverB_C4Grass", "prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon") %>% drop_na(),
response_vars = c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), # name of observations, followed by name of predictions
pred_vars = c("prcp", "prcp_dry", "prcpTempCorr", "isothermality", "sand", "coarse", "AWHC", "isothermality_anom", "tmean", "prcp_anom", "clay", "carbon"), # for predictors, combine list of predictors for all models (total tree for GS and F, and perc BL tree for GS and F)
cut_points = seq(0, 1, 0.01)
)# shrubs
quantPlot_shrubs_Obs <- decile_dotplot_pq(df = deciles_shrub_Obs, response= c("relCoverB_shrub_Obs", "relCoverB_shrub"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Shrub Cover: Quantile Plot")
quantPlot_shrubs_Obs <- add_dotplot_inset(quantPlot_shrubs_Obs, df = deciles_shrub_Obs, response = c("relCoverB_shrub_Obs", "relCoverB_shrub"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# bare ground
quantPlot_bareGround_Obs <- decile_dotplot_pq(df = deciles_bareGround_Obs, response= c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Bare Ground Cover: Quantile Plot")
quantPlot_bareGround_Obs <- add_dotplot_inset(quantPlot_bareGround_Obs, df = deciles_bareGround_Obs, response = c("relCoverB_bareGround_Obs", "relCoverB_bareGround"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# total herbaceous
quantPlot_totHerb_Obs <- decile_dotplot_pq(df = deciles_totalHerb_Obs, response= c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Total Herbaceous Cover: Quantile Plot")
quantPlot_totHerb_Obs <- add_dotplot_inset(quantPlot_totHerb_Obs, df = deciles_totalHerb_Obs, response = c("relCoverB_totalHerb_Obs", "relCoverB_totalHerb"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# forbs
quantPlot_forbs_Obs <- decile_dotplot_pq(df = deciles_forb_Obs, response= c("relCoverB_Forb_Obs", "relCoverB_Forb"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute Forb Cover: Quantile Plot")
quantPlot_forbs_Obs <- add_dotplot_inset(quantPlot_forbs_Obs, df = deciles_forb_Obs, response = c("relCoverB_Forb_Obs", "relCoverB_Forb"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
# C3 grass
quantPlot_c3grass_Obs <- decile_dotplot_pq(df = deciles_C3grass_Obs, response= c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute C3 Grass Cover: Quantile Plot")
quantPlot_c3grass_Obs <- add_dotplot_inset(quantPlot_c3grass_Obs, df = deciles_C3grass_Obs, response = c("relCoverB_C3Grass_Obs", "relCoverB_C3Grass"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)
#C4 grass
quantPlot_c4grass_Obs <- decile_dotplot_pq(df = deciles_C4grass_Obs, response= c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), IQR = TRUE,
CI = FALSE
) + ggtitle("Absolute C4 Grass Cover: Quantile Plot")
quantPlot_c4grass_Obs <- add_dotplot_inset(quantPlot_c4grass_Obs, df = deciles_C4grass_Obs, response = c("relCoverB_C4Grass_Obs", "relCoverB_C4Grass"), dfRaw = preds_quantile, add_smooth = TRUE, deciles = FALSE)## now, display the figures
quantPlot_totHerb_Obs
quantPlot_shrubs_Obs
quantPlot_bareGround_Obs
quantPlot_c3grass_Obs
quantPlot_c4grass_Obs
quantPlot_forbs_Obs
Below, we show the RMSE and bias of the final estimates of relative cover for each functional type (except trees) in comparison to the observed data.
long_absPreds <- preds_quantile %>%
mutate(uniqueID = 1:nrow(preds_quantile)) %>%
select(uniqueID, relCoverB_totalHerb, relCoverB_shrub,
relCoverB_bareGround, relCoverB_C3Grass,
relCoverB_C4Grass, relCoverB_Forb) %>%
rename(TotalHerb = relCoverB_totalHerb,
Shrub = relCoverB_shrub,
BareGround = relCoverB_bareGround,
C3gram = relCoverB_C3Grass,
C4gram = relCoverB_C4Grass,
Forb = relCoverB_Forb) %>%
pivot_longer(cols = c(TotalHerb, Shrub, BareGround, C3gram, C4gram, Forb),
names_to = "names",
values_to = "preds_values"
) %>%
select(uniqueID, names, preds_values)
long_absObs <- preds_quantile %>%
mutate(uniqueID = 1:nrow(preds_quantile)) %>%
select(uniqueID, relCoverB_totalHerb_Obs, relCoverB_shrub_Obs, relCoverB_bareGround_Obs, relCoverB_C3Grass_Obs, relCoverB_C4Grass_Obs, relCoverB_Forb_Obs) %>%
rename(TotalHerb = relCoverB_totalHerb_Obs,
Shrub = relCoverB_shrub_Obs,
BareGround = relCoverB_bareGround_Obs,
C3gram = relCoverB_C3Grass_Obs,
C4gram = relCoverB_C4Grass_Obs,
Forb = relCoverB_Forb_Obs) %>%
pivot_longer(cols = c(TotalHerb, Shrub, BareGround, C3gram, C4gram, Forb),
names_to = "names",
values_to = "obs_values"
) %>%
select(uniqueID, names, obs_values)
longObsPreds <- long_absObs %>%
left_join(long_absPreds)
## calculate bias
(summaryTableAbs <- longObsPreds %>%
mutate(resids = obs_values - preds_values) %>%
group_by(names) %>%
dplyr::summarise(bias = mean(resids, na.rm = TRUE),
RMSE = round(sqrt(mean(resids^2, na.rm = TRUE)),4)
) %>%
slice(8, 9, 7, 2, 6, 1, 3, 4, 5) %>% kable(
format = "html",
col.names = c("Model Name", "Model bias mean(obs. - pred.)", "model RMSE")
))| Model Name | Model bias mean(obs. - pred.) | model RMSE |
|---|---|---|
| C3gram | 0.0207554 | 0.1798 |
| TotalHerb | -0.0136944 | 0.2156 |
| BareGround | 0.0219762 | 0.1887 |
| C4gram | 0.0012859 | 0.1014 |
| Forb | -0.0146440 | 0.1428 |
| Shrub | -0.0082818 | 0.1760 |
Finally, we show how these model predictions made using contemporary climate and weather data change over time from 2010 to 2020.
# calculate relativized cover (version w/ all functional groups summing to 1)
predsOverTime <- predsOverTime %>%
mutate(total_NonTree_AbsCover_CONUS = (absTotalHerb_CONUS + #absTotalTree_CONUS +
absBareGround_CONUS + absShrub_CONUS)) %>%
mutate(relCoverB_totalTree = absTotalTree_CONUS,#/total_NonTree_AbsCover_CONUS,
relCoverB_totalHerb = absTotalHerb_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_shrub = absShrub_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_bareGround = absBareGround_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_NLTree = absNLTree_CONUS,
relCoverB_BLTree = absBLTree_CONUS,
relCoverB_C3Grass = absC3grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_C4Grass = absC4grass_CONUS/total_NonTree_AbsCover_CONUS,
relCoverB_Forb = absForb_CONUS/total_NonTree_AbsCover_CONUS
)
# make into a long format
predsOverTime_long <- predsOverTime %>%
select(year,x, y, spatialID, uniqueID, total_NonTree_AbsCover_CONUS:relCoverB_Forb) %>%
pivot_longer(cols = c(total_NonTree_AbsCover_CONUS:relCoverB_Forb)) %>%
mutate(year = as.integer(year))
# calculate the coefficient of variation for each cover type at each location
predsOverTime_cv <- predsOverTime_long %>%
group_by(uniqueID, name) %>%
summarize(cv = sd(value)/mean(value))
## add x,y
predsOverTime_cv <- predsOverTime_cv %>%
left_join(predsOverTime %>% select(uniqueID, x, y))
# visualize change over time for a selection of points
coverOverTimeForSelectLocations <- predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c( #"total_NonTree_AbsCover_CONUS", "relCoverB_totalTree",
#"relCoverB_totalHerb",
"relCoverB_shrub" , "relCoverB_bareGround" , "relCoverB_NLTree", "relCoverB_BLTree",
"relCoverB_C3Grass" , "relCoverB_C4Grass", "relCoverB_Forb" )) %>%
ggplot() +
facet_wrap(~uniqueID) +
geom_area(aes(x = year, y = value, fill = name)) +
geom_line(data = predsOverTime_long %>%
filter(uniqueID %in% 1:50) %>%
filter(name %in% c("relCoverB_totalHerb", "relCoverB_totalTree")),
aes(x = year, y = value, linetype = name), col = "black") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90)) +
labs(title = "Change in modeled relative cover (all Non-tree types sum to 1) by functional type from \n 2010 to 2020 for a subset of locations")
mapOfSelectLocations <- predsOverTime_long %>%
select(x, y, uniqueID) %>%
filter(uniqueID %in% 1:50) %>%
unique() %>%
ggplot() +
geom_sf(data = cropped_states, aes()) +
geom_point(aes(x, y)) +
geom_label_repel(aes(x, y, label = uniqueID)) +
labs(title = "Locations in above plot")
ggarrange(coverOverTimeForSelectLocations, mapOfSelectLocations, ncol = 1, heights = c(1,.5))
predsOverTime_cv %>%
ggplot() +
geom_boxplot(aes(x = name, y = cv, col = name)) +
theme_minimal() +
labs(title = "Coefficients of variation of model-predicted relative cover (all Non-tree types sum to 1) \n at a given location for each functional group from 2010-2020",
x = "functional group",
y = "coefficient of variation") +
theme(axis.text.x = element_text(angle = 90))
# show the locations of each point that is being used for this temporal analysis
#cropped_states <- cropped_states %>%
# rename(geometry = Shape)
ggplot() +
facet_wrap(~name) +
geom_sf(data = cropped_states, aes()) +
stat_summary2d(data = predsOverTime_cv, aes(x = x, y = y, z = cv), fun = ~ mean(.x, na.rm = TRUE)) +
geom_hex(data = predsOverTime_cv, aes(x = x, y = y, col = cv), stat = "mean", na.rm = TRUE) +
guides(fill = guide_legend(title = "CV")) +
theme_minimal() +
labs(title = "The spatial distribution of the coefficients of variation of \nmodel-predicted cover from 2010-2020 for 1000 random locations",
subtitle = "Note that values are aggregated accross space to enhance readability")